#include "stdafx.h"
#include "CroppedTask.h"

CroppedTask::CroppedTask(int x, int y, double acc, int steps): Task(x,y,acc,steps)
{
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
	{
		matrix[i]=0.0;
		matrix2[i]=0.0;
		diff[i]=0.0;
	}
}

void CroppedTask::solve()
{
	double*	f = new double[(size_x + 1)*(size_y + 1)];
	double x,y;
	bool flag = false;
	double eps_max = 0;
	double eps_cur = 0;
	int steps = 0;

	int m_x = (size_x + 1)/2;
	int m_y = (size_y + 1)/2;

	double h2 = size_x*size_x;
	double k2 = size_y*size_y;
	double a2 = -2 * (h2 + k2);

	for(int j=0; j<size_y+1; j++)
		for (int i=0; i<size_x+1; i++)
		{
			x = 1.0 + (1.0/size_x)*i; 
			y = 1.0 + (1.0/size_y)*j;
			f[(size_x + 1)*j + i] = 4*x*x*exp(x*x - y*y) + 4*y*y*exp(x*x - y*y);
		}

	// up
	for (int i=0; i<size_x+1; i++)
	{
		x = 1.0 + 1.0*i/size_x; 
		matrix[size_y*(size_x + 1) + i] = exp(x*x - 4.0);
		f[(size_y - 1)*(size_x + 1) + i] -= h2*exp(x*x - 4.0); 
	}
	// down1
	for (int i=0; i<(size_x+1)/2; i++)
	{
		x = 1.0 + 1.0*i/size_x; 
		matrix[i] = exp(x*x - 1.0);
		f[i+size_x+1] -= h2*exp(x*x - 1.0);
	}
	//down2
	for (int i=(size_x+1)/2; i<size_x+1; i++)
	{
		x = 1.0 + 1.0*i/size_x;
		y = 1.0 + 1.0*((double) m_y)/size_y;
		matrix[i+((size_y+1)/2)*(size_x+1)] = exp(x*x - y*y);
		f[i+((size_y+1)/2 + 1)*(size_x+1)] -= h2*exp(x*x - y*y);
	}
	//left
	for (int i=0; i<size_y+1; i++)
	{
		y = 1.0 + 1.0*i/size_y;
		matrix[(size_x + 1)*i] = exp(1 - y*y);
		f[(size_x + 1)*i + 1] -= k2*exp(1 - y*y);
	}
	//right1
	for (int i=0; i<(size_y+1)/2+1; i++)
	{
		x = 1.0 + 1.0*((double) m_x)/size_x;
		y = 1.0 + 1.0*i/size_y;
		matrix[(size_x + 1)*i + (size_y+1)/2] = exp(x*x - y*y);
		f[(size_x + 1)*i + (size_y+1)/2 - 1] -= k2*exp(x*x - y*y);
	}
	//right2
	for (int i=(size_y+1)/2; i<size_y+1; i++)
	{
		y = 1.0 + 1.0*i/size_y;
		matrix[(size_x + 1)*i + size_y] = exp(4 - y*y);
		f[(size_x + 1)*i + size_y - 1] -= k2*exp(4 - y*y);
	}
		
	// start method
	double* x_vector = new double[(size_x + 1)*(size_y + 1)];
	double* h = new double[(size_x + 1)*(size_y + 1)];
	double* r = new double[(size_x + 1)*(size_y + 1)];
	double beta = 0.0;
	double alpha = 0.0;

	// Initialize matrix of zeros
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
	{
		r[i] = 0.0;
		x_vector[i] = 0.0;
		h[i] = 0.0;
	}
	r = updateR(x_vector, r, f);
	// h = -r;
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
		h[i] = -r[i];
	alpha = getAlpha(h,r);

	while (!flag)
	{
		eps_max = 0.0;
		// update x
		for (int j=1; j<(size_y+1)/2+1; j++ )
			for (int i=1; i<(size_x+1)/2; i++)
			{
				x_vector[(size_x + 1)*j + i] = x_vector[(size_x + 1)*j + i] - alpha*h[(size_x + 1)*j + i];
				eps_cur = fabs(alpha*h[(size_x + 1)*j + i]);
				if (eps_cur > eps_max)
					eps_max = eps_cur;
			}
		for (int j=(size_y+1)/2+1; j<size_y; j++ )
			for (int i=1; i<size_x; i++)
			{
				x_vector[(size_x + 1)*j + i] = x_vector[(size_x + 1)*j + i] - alpha*h[(size_x + 1)*j + i];
				eps_cur = fabs(alpha*h[(size_x + 1)*j + i]);
				if (eps_cur > eps_max)
					eps_max = eps_cur;
			}
		steps++;
		if ((eps_max < accuracy) || (steps >= num_steps)) flag = true;
		// calc r(s)
		r = updateR(x_vector, r, f);
		beta = getBeta(h,r);
		if (fabs(beta) < beta_accuracy) {
			flag = true;
			break;
		}
		// calc h(s)
		h = getH(h, r, beta);
		//calc alpha(s)
		alpha = getAlpha(h,r);
	}
	accuracy = eps_max;
	num_steps = steps;

	// store result
	// down
	for (int j=1; j<(size_y+1)/2+1; j++ )
		for (int i=1; i<(size_x+1)/2; i++)
			matrix[(size_x + 1)*j + i] = x_vector[(size_x + 1)*j + i];
	// up
	for (int j=(size_y+1)/2+1; j<size_y; j++ )
		for (int i=1; i<size_x; i++)
			matrix[(size_x + 1)*j + i] = x_vector[(size_x + 1)*j + i];

	// calculate matrix2 and difference
	// down
	for (int j=0; j<(size_y+1)/2; j++ )
		for (int i=0; i<(size_x+1)/2+1; i++)
		{
			x = 1.0 + 1.0/size_x*i; 
			y = 1.0 + 1.0/size_y*j;
			matrix2[(size_x + 1)*j + i] = exp(x*x - y*y);
			diff[(size_x + 1)*j + i] = fabs(matrix[(size_x + 1)*j + i] - matrix2[(size_x + 1)*j + i]);
		}
	// up
	for (int j=(size_y+1)/2; j<size_y+1; j++ )
		for (int i=0; i<size_x+1; i++)
		{
			x = 1.0 + 1.0/size_x*i; 
			y = 1.0 + 1.0/size_y*j;
			matrix2[(size_x + 1)*j + i] = exp(x*x - y*y);
			diff[(size_x + 1)*j + i] = fabs(matrix[(size_x + 1)*j + i] - matrix2[(size_x + 1)*j + i]);
		}

	// release arrays
	delete[] h;
	delete[] r;
	delete[] x_vector;
}

double* CroppedTask::updateR(double* x, double* r, double* f)
{
	double* temp = new double[(size_x + 1)*(size_y + 1)];
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
		temp[i] = x[i];

	temp = mulA(temp);
	for (int j=1; j<(size_y+1)/2+1; j++ )
		for (int i=1; i<(size_x+1)/2; i++)
			r[(size_x + 1)*j + i] = temp[(size_x + 1)*j + i] - f[(size_x + 1)*j + i];
	
	for (int j=(size_y+1)/2+1; j<size_y; j++ )
		for (int i=1; i<size_x; i++)
			r[(size_x + 1)*j + i] = temp[(size_x + 1)*j + i] - f[(size_x + 1)*j + i];

	delete[] temp;
	return r;
}

double CroppedTask::getBeta(double* h, double* r)
{
	//input: h - h(s-1), r - r(s) 
	double beta = 0.0;
	// get copy h
	double* temp = new double[(size_x + 1)*(size_y + 1)];
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
		temp[i] = h[i];

	temp = mulA(temp);
	beta = scalarMul(temp, r) / scalarMul(temp, h);
	delete[] temp;
	return beta;
}

double* CroppedTask::getH(double* h, double* r, double beta)
{
	// input: h - h(s-1), r - r(s)  
	//double beta = getBeta(h,r);
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
		h[i] = -r[i] + beta*h[i];
	return h;
}

double CroppedTask::getAlpha(double* h, double* r)
{
	// input: h - h(s-1), r - r(s) 
	double alpha = 0.0;
	// get copy h
	double* temp = new double[(size_x + 1)*(size_y + 1)];
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
		temp[i] = h[i];

	temp = mulA(temp);
	alpha = scalarMul(r,h) / scalarMul(temp, h);
	delete[] temp;
	return alpha;
}

double* CroppedTask::mulA(double* matrix)
{
	double* temp = new double[(size_x + 1)*(size_y + 1)];
	for (int i=0; i<(size_x+1)*(size_y+1); i++)
		temp[i] = matrix[i];
	double h2 = size_x*size_x;
	double k2 = size_y*size_y;
	double a2 = -2 * (h2 + k2);
	// down
	for (int j=1; j<(size_y+1)/2+1; j++ )
		for (int i=1; i<(size_x+1)/2; i++)
			matrix[(size_x + 1)*j + i] = k2*(temp[(size_x + 1)*(j+1) + i] + temp[(size_x + 1)*(j-1) + i]) + h2*(temp[(size_x + 1)*j + i+1] + temp[(size_x + 1)*j + i-1]) + temp[(size_x + 1)*j + i]*a2;
	// up
	for (int j=(size_y+1)/2+1; j<size_y; j++ )
		for (int i=1; i<size_x; i++)
			matrix[(size_x + 1)*j + i] = k2*(temp[(size_x + 1)*(j+1) + i] + temp[(size_x + 1)*(j-1) + i]) + h2*(temp[(size_x + 1)*j + i+1] + temp[(size_x + 1)*j + i-1]) + temp[(size_x + 1)*j + i]*a2;
	delete[] temp;
	return matrix;
}

double CroppedTask::scalarMul(double* v1, double* v2)
{
	double result = 0.0;
	// down
	for (int j=1; j<(size_y+1)/2+1; j++ )
		for (int i=1; i<(size_x+1)/2; i++)
			result+=v1[(size_x + 1)*j + i]*v2[(size_x + 1)*j + i];
	// up
	for (int j=(size_y+1)/2+1; j<size_y; j++ )
		for (int i=1; i<size_x; i++)
			result+=v1[(size_x + 1)*j + i]*v2[(size_x + 1)*j + i];
	return result;
}